## 

library(tidyverse)
library(lfe)
library(pbapply)
library(rdrobust)
library(fixest)
library(broom)

## ##

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')
drop_state <- '08'

## Get data - municipal elections 

mu <- readRDS('data/data_municipal.rds') %>%
  filter(!is.na(treated)) %>%
  mutate(state_id = substr(ags, 1, 2)) %>% 
  filter(time_rel_period > -2) %>% 
  filter(!state_id == '13') %>% ## This state has just one election in the data, so discard
  filter(!state_id == '08')   ### BW excluded

## Check election dates by state 

state_id_to_names <- function (state_id) 
{
  names <- recode(state_id, `01` = "Schleswig-Holstein", `02` = "Hamburg", 
                  `03` = "Niedersachsen", `04` = "Bremen", `05` = "North Rhine-Westphalia", 
                  `06` = "Hesse", `07` = "Rhineland-Palatinate", `08` = "Baden-Württemberg", 
                  `09` = "Bavaria", `10` = "Saarland", `11` = "Berlin", 
                  `12` = "Brandenburg", `13` = "Mecklenburg-Vorpommern", 
                  `14` = "Saxony", `15` = "Saxony-Anhalt", `16` = "Thuringia")
  return(names)
}

##

mu %>% group_by(state_id) %>% 
  summarise(ed = paste0(unique(year) %>% sort(), collapse = ', ')) %>% 
  mutate(state_id = state_id_to_names(state_id))

## Get unem data

unem <- read_rds("data/muni_unem.rds") %>% 
  dplyr::select(ags, year, unem_capita)

## Get county data

unem_county <- read_rds("data/county_wage_gdp_unem.rds") %>% 
  dplyr::select(ags, year, gdp_nominal_capita, unem_rate_tot) %>% 
  dplyr::rename(county_gdp_nominal_capita = gdp_nominal_capita,
                county_unem_rate_tot = unem_rate_tot) %>% 
  dplyr::rename(ags_county = ags)

## Merge this to the main data

mu <- mu %>% 
  left_join(unem) %>% 
  left_join(unem_county)

## Code change in unem

mu <- mu %>% 
  arrange(ags, year) %>% 
  group_by(ags) %>% 
  mutate(unem_rate_change = (unem_capita[2] - unem_capita[1])*100,
         unem_rate_increase = ifelse(unem_rate_change > 0, 1, 0),
         county_unem_rate_change = (county_unem_rate_tot[2] - 
                                      county_unem_rate_tot[1])*100,
         county_unem_rate_increase = ifelse(county_unem_rate_change > 0, 1, 0),
         county_gdp_change = (county_gdp_nominal_capita[2] - 
                                county_gdp_nominal_capita[1])/1000,
         county_gdp_increase = ifelse(county_gdp_change > 0, 1, 0))

## Relevant vars

# inc_share2 - incumbent share (only for the period right before and after census)
# inc_party_prior_to_census - incumbent prior to census

## Outcomes

outcomes <- c('inc_share2')

## Get differenced DF

olist <- c("inc_share2", "agg_left")

diff_df <- pblapply(olist, function(o) {
  out <- mu %>%
    filter(time_rel_period > -2) %>%
    mutate(time_rel_period = time_rel_period + 2)  %>% 
    pivot_wider(values_from = all_of(o), names_from = 'time_rel_period', 
                id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2  - o1) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(mu %>% dplyr::select(ags, pop_dec_09, 
                                 matches('unem_rate|county'),
                                 one_of(covs), state_id,
                                 inc_party_prior_to_census,
                                 applies_census) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1)

## Estimate main model to get BW

m_main <- rdrobust(y = diff_df[, 'agg_left'] %>% pull(agg_left), 
                   x = diff_df$runvar, c = 0,
                   covs = diff_df[, covs])
bw_use <- m_main$bws[1, 1]

## Regressions
## Muni level

library(fixest)

m1 <- feols(inc_share2 ~ I(scale(unem_rate_change)), data = diff_df %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(inc_share2 ~ I(scale(unem_rate_change)) | state_id, data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(inc_share2 ~ unem_rate_increase, data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(inc_share2 ~ unem_rate_increase | state_id, data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

out_muni_munitreat <- list(m1, m2, m3, m4) %>% 
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'unem_rate')) %>% 
  mutate(unit = 'Incumbent in\nmunicipal elections', treat_unit = 'muni',
         fe = rep(c("none", "state"), 2))

## County-level regressions ##

m1 <- feols(inc_share2 ~ I(scale(county_unem_rate_change)), data = diff_df %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(inc_share2 ~ I(scale(county_unem_rate_change)) | state_id, 
            data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(inc_share2 ~ I(scale(county_gdp_change)), data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(inc_share2 ~ I(scale(county_gdp_change)) | state_id, 
            data = diff_df%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

out_muni_countytreat <- list(m1, m2, m3, m4)%>% 
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'county')) %>% 
  mutate(unit = 'Incumbent in\nmunicipal elections', treat_unit = 'county',
         fe = rep(c("none", "state"), 2))

## Federal elections ##

bt <- readRDS('data/data_federal.rds') %>%
  filter(!is.na(treated)) %>% 
  filter(year > 2012) %>% 
  dplyr::select(-unem_capita)

## Prep

olist <- c("incumbent_party")

## Merge this to the main data

bt <- bt %>% 
  left_join(unem) %>% 
  left_join(unem_county)

## Code change in unemployment

bt <- bt %>% 
  arrange(ags, year) %>% 
  group_by(ags) %>% 
  mutate(unem_rate_change = (unem_capita[2] - unem_capita[1])*100,
         unem_rate_increase = ifelse(unem_rate_change > 0, 1, 0),
         county_unem_rate_change = (county_unem_rate_tot[2] - 
                                      county_unem_rate_tot[1])*100,
         county_unem_rate_increase = ifelse(county_unem_rate_change > 0, 1, 0),
         county_gdp_change = (county_gdp_nominal_capita[2] - 
                                county_gdp_nominal_capita[1])/1000,
         county_gdp_increase = ifelse(county_gdp_change > 0, 1, 0))

## Relevant vars

# incumbent_party - incumbent share (only for the period right before and after census)
# inc_party_prior_to_census - incumbent prior to census

## First differences

olist <- c("incumbent_party", "agg_left_party")

diff_df_bt <- pblapply(olist, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id,
                                 matches("unem|gdp")) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1) %>% 
  filter(!state_id == '08')

## Estimate main model to get BW

m_main <- rdrobust(y = diff_df_bt[, 'agg_left_party'] %>% pull(agg_left_party), 
                   x = diff_df_bt$runvar, c = 0,
                   covs = diff_df_bt[, covs])
bw_use <- m_main$bws[1, 1]

## Regressions
## Muni level

m1 <- feols(incumbent_party ~ I(scale(unem_rate_change)), data = diff_df_bt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(incumbent_party ~ I(scale(unem_rate_change)) | state_id, data = diff_df_bt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(incumbent_party ~ unem_rate_increase, data = diff_df_bt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(incumbent_party ~ unem_rate_increase | state_id, data = diff_df_bt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

## ## 

out_bt_munitreat <- list(m1, m2, m3, m4) %>%  
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'unem_rate')) %>% 
  mutate(unit = 'Incumbent in\nfederal elections',
         treat_unit = 'muni',
         fe = rep(c("none", "state"), 2))

## County-level regressions ##

m1 <- feols(incumbent_party ~ I(scale(county_unem_rate_change)), data = diff_df_bt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(incumbent_party ~ I(scale(county_unem_rate_change)) | state_id, 
            data = diff_df_bt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(incumbent_party ~ I(scale(county_gdp_change)), data = diff_df_bt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(incumbent_party ~ I(scale(county_gdp_change)) | state_id, 
            data = diff_df_bt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

out_bt_countytreat <- list(m1, m2, m3, m4) %>%  
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'unem_rate|gdp'))  %>% 
  mutate(unit = 'Incumbent in\nfederal elections', 
         treat_unit = 'county',
         fe = rep(c("none", "state"), 2))

## State elections ##

lt <- readRDS('data/data_state.rds') %>%
  filter(!is.na(treated)) %>%
  mutate(year = lubridate::year(date)) %>%
  mutate(state_id = substr(ags, 1, 2)) %>% 
  dplyr::select(-unem_capita)

## Get list of states to use

lt <- read_rds("data/states_census.rds") %>%
  dplyr::select(applies_census, state_id, census_first_year) %>% 
  left_join(lt, .) %>%
  mutate(time_rel = year - census_first_year)

## This converts each election year to time period relative to treatment
## Errors can be ignored

periods_by_state <- lt %>%
  group_by(state_id) %>%
  summarise(l = list(unique(time_rel))) %>%
  unnest()
periods_by_state <- periods_by_state %>%
  group_by(state_id) %>%
  summarise(l_r = list(rank(l))) %>%
  unnest() %>% dplyr::select(-state_id) %>% 
  bind_cols(periods_by_state, .)
periods_by_state <- periods_by_state %>% 
  group_by(state_id) %>%
  mutate(l_r = l_r - max(l_r)) %>%
  filter(!is.na(l)) %>%
  rename(time_rel_period = l_r,
         time_rel = l)

## Merge this back to the main DF

lt <- lt %>%
  left_join(periods_by_state) %>%
  mutate(post = ifelse(time_rel_period > -1, 1, 0),
         treated_post = post * treated) %>% 
  filter(time_rel_period > -2)

## Merge unem to the main data

lt <- lt %>% 
  left_join(unem) %>% 
  left_join(unem_county)

## Code change in unem

lt <- lt %>% 
  arrange(ags, year) %>% 
  group_by(ags) %>% 
  mutate(unem_rate_change = (unem_capita[2] - unem_capita[1])*100,
         unem_rate_increase = ifelse(unem_rate_change > 0, 1, 0),
         county_unem_rate_change = (county_unem_rate_tot[2] - 
                                      county_unem_rate_tot[1])*100,
         county_unem_rate_increase = ifelse(county_unem_rate_change > 0, 1, 0),
         county_gdp_change = (county_gdp_nominal_capita[2] - 
                                county_gdp_nominal_capita[1])/1000,
         county_gdp_increase = ifelse(county_gdp_change > 0, 1, 0))

## First differences 

olist <- c("agg_left", "incumbent")

diff_df_lt <- pblapply(olist, function(o) {
  out <- lt %>%
    filter(time_rel_period > -2) %>%
    mutate(time_rel_period = time_rel_period + 2)  %>% 
    pivot_wider(values_from = o, names_from = 'time_rel_period', 
                id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2  - o1) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(lt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs),
                                 state_id,
                                 matches("unem|county")) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1) %>%
  filter(!str_sub(ags, 1, 2) == '01') %>% ## Exclude SH
  filter(!(substr(ags, 1, 2) == '08')) ## Exclude BW

## Estimate main model to get BW

m_main <- rdrobust(y = diff_df_lt[, 'agg_left'] %>% pull(agg_left), 
                   x = diff_df_lt$runvar, c = 0,
                   covs = diff_df_lt[, covs])
bw_use <- m_main$bws[1, 1]

## Regressions
## Muni level

m1 <- feols(incumbent ~ I(scale(unem_rate_change)), data = diff_df_lt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(incumbent ~ I(scale(unem_rate_change)) | state_id, data = diff_df_lt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(incumbent ~ unem_rate_increase, data = diff_df_lt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(incumbent ~ unem_rate_increase | state_id, data = diff_df_lt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

## ## 

out_lt_munitreat <- list(m1, m2, m3, m4) %>%  
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'unem_rate')) %>% 
  mutate(unit = 'Incumbent in\nstate elections', treat_unit = 'muni',
         fe = rep(c("none", "state"), 2))

## County-level regressions ##

m1 <- feols(incumbent ~ I(scale(county_unem_rate_change)), data = diff_df_lt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m2 <- feols(incumbent ~ I(scale(county_unem_rate_change)) | state_id, 
            data = diff_df_lt%>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m3 <- feols(incumbent ~ I(scale(county_gdp_change)), data = diff_df_lt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))
m4 <- feols(incumbent ~ I(scale(county_gdp_change)) | state_id, 
            data = diff_df_lt %>% 
              filter(between(pop_dec_09, 10000-bw_use, 10000+bw_use)))

## ## 

out_lt_countytreat <- list(m1, m2, m3, m4) %>%  
  lapply(broom::tidy) %>% 
  reduce(rbind) %>% 
  filter(str_detect(term, 'unem_rate|gdp')) %>% 
  mutate(unit = 'Incumbent in\nstate elections', treat_unit = 'county',
         fe = rep(c("none", "state"), 2))

rm(list = ls()[!str_detect(ls(), "out_|diff_df")])

## Combine results

out_full <- list(out_bt_countytreat, 
                 out_bt_munitreat, 
                 out_lt_countytreat, 
                 out_lt_munitreat, 
                 out_muni_countytreat, 
                 out_muni_munitreat) %>% 
  reduce(rbind)

## Dict for treatments

dict_df <- data.frame(term = c("I(scale(county_unem_rate_change))", "I(scale(county_gdp_change))", 
                               "I(scale(unem_rate_change))", "unem_rate_increase"),
                      treat = c("Change in unemployment rate\n(county-level, standardized)",
                                "Change in GDP/capita\n(county-level, standardized)",
                                "Change in unemployment rate\n(municipality-level, standardized)",
                                "Any increase in unemployment rate\n(municipality-level, dichotomous)"),
                      order = c(3,4,1,2))

dict_df_unit <- data.frame(unit = unique(out_full$unit),
                           order_unit = c(3,2,1))

## merge

out_full <- out_full %>% 
  left_join(dict_df) %>% 
  left_join(dict_df_unit) %>% 
  mutate(treat = factor(treat),
         treat = fct_reorder(treat, order, .desc = T)) %>% 
  mutate(unit = factor(unit),
         unit = fct_reorder(unit, order_unit, .desc = T)) 

# Figure A.25: economic conditions and incumbent support ----

out_full %>%   filter(fe == 'state') %>% 
  ggplot(aes(treat, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, 
                    ymax = estimate + 1.96*std.error),
                width  = 0) +
  geom_point(shape = 21, fill = 'white', size = 2) +
  facet_wrap(~unit) +
  coord_flip() +
  theme_bw() +
  ylab("Change in incumbent vote share\n(percentage points)") +
  xlab("Independent variable")

